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1.0 


INTRODUCTION 


Nonlinear  structural  analyses  include  materially  nonlinear  problems  such  as  nonlinear  constitutive 
equations  with  small  deformations,  geometrically  nonlinear  problems  normally  associated  with 
buckling,  or  a  combination  of  both  types  of  nonlinearities  (1).  If  the  tangent  stiffness  matrix,  K,  is 
defined  to  include  both  types  of  nonlinearities,  a  limit  or  bifurcation  point  is  characterized  by  a 
critical  state  at  which  the  tangent  stiffness  matrix  satisfies  the  static  stability  criterion 

KcAuc  =  0  (1-1) 

in  which  Auc  denotes  the  vector  (eigenmode)  of  displacement  increments.  Solutions  for  the 
displacement  vector  beyond  critical  points  are  important  for  describing  the  failure  process  and  for 
making  comparisons  with  experimental  data  (2, 3). 

Most  solution  schemes  for  nonlinear  problems  are  based  on  step-by-step  load  incrementation  and 
an  iteration  procedure  to  correct  for  the  linearization  error.  To  trace  nonlinear  response  beyond 
critical  points,  several  methods  have  been  proposed  to  circumvent  the  singularity  that  appears  in  the 
tangent  stiffness  matrix.  One  method  is  to  detect  the  presence  of  negative  pivots  from  the 

decomposition  of  Kc  and  then  replace  Kc  with  an  expression  of  the  type  Kc  +  rjl  with  t|  >  0  so 
that  this  matrix  is  "safely"  positive  definite  in  the  vicinity  of  critical  points  (4).  A  second  method 
consists  of  suppressing  the  iterations  around  the  critical  point  and  reversing  the  sign  of  the  load 
following  the  appearance  of  the  negative  pivot  (5, 6).  In  order  to  avoid  deviating  too  far  from  the 
equilibrium  path,  very  small  load  steps  have  to  be  used  in  the  range  where  no  iterations  are 
performed.  Because  it  is  extremely  difficult  to  obtain  converged  solutions  around  critical  points  by 
simply  prescribing  external  load  increments,  it  has  been  suggested  that  the  load  level  become  a 
function  of  another  variable.  Direct  displacement  control  (7),  arc-length  control  (8),  self-adaptive 
"hyperelastic"  constraints  (9),  and  indirect  displacement  control  (10)  are  designed  to  provide  such  a 
constraint  function  through  which  the  load  level  changes  from  iteration  to  iteration.  For  the  case  of 
snap-through  in  geometrically  nonlinear  problems,  direct  displacement  control  fails  and  arc-length 
control  appears  to  be  a  robust  procedure.  However,  for  materially  nonlinear  analyses,  a  global 
norm  on  incremental  displacements  as  used  in  arc -control  is  often  less  successful  due  to  localization 
effects,  and  it  may  be  more  appropriate  to  employ  only  one  dominant  degree  of  freedom  or  to  omit 
some  degrees  of  freedom  from  the  norm  of  incremental  displacements.  Such  a  scheme  is  often 
referred  to  as  one  of  "indirect  displacement  control."  The  disadvantage  of  modifying  the  constraint 
condition  is  that  the  constraint  function  becomes  problem  dependent. 
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Recently,  with  the  increase  in  interest  in  softening  which  is  accompanied  by  localization,  the  need 
for  a  robust  solution  strategy  has  become  of  paramount  importance.  Efforts  have  been  made 
towards  modifying  arc-length  control  so  that  the  constraint  equation  is  sensitive  to  the  change  of 
state  variables  w  ithin  the  localization  zone.  One  modification  is  to  use  the  eigenvector  of  the  lowest 
eigenvalue  of  the  tangent  stiffness  matrix  as  a  trigger  to  move  the  solution  from  the  unstable, 
unlocalized  path  (11, 12).  Another  one  is  the  minimization  of  a  function  of  an  out-of-balance  force 
vector  to  seek  a  dominant  direction  and  a  line  search  procedure  along  that  direction  (4).  The  line 
search  method  improves  the  convergence  characteristics  of  arc-length  control  by  introducing  a 
variable  step-length  parameter,  which  scales  the  usual  iterative  displacement  vector  (13).  Since 
these  solution  strategies  have  not  been  applied  routinely  to  multidimensional  problems  involving 
general  constitutive  equations,  they  must  be  considered  as  tentative  measures  in  the  search  for  a 
general  robust  scheme.  A  better  mathematical  foundation  is  required  for  the  application  of  line 
search  techniques  when  tracing  unstable  equilibrium  paths,  and  the  issue  of  appropriate  incremental 
quantities  needs  further  consideration.  It  has  been  pointed  out  that  existing  procedures  will  not 
eliminate  the  difficulties  associated  with  the  dependence  on  increment  sizes  and  on  alternative 
equilibrium  states,  although  nonlocal  constitutive  models  may  reduce  the  mesh-dependency  of 
numerical  solutions  (11). 

In  general,  arc-length  control  is  still  commonly  used  in  geometrically  nonlinear  problems  and,  with 
some  modifications,  in  materially  nonlinear  analyses  since  arc-length  control  can  trace  the  post-peak 
response  even  for  snap-through  and  snap-back  problems.  Occasionally,  lack  of  convergence  is 
reported  even  after  the  limit  or  bifurcation  point  has  been  traversed.  One  reason  might  be  due  to  the 
fact  that  the  load  increment  in  the  first  iteration  has  to  be  prescribed  according  to  some  condition, 
and  there  is  a  sign  change  of  the  load  increment  during  the  transition  from  the  pre-peak  to  the  post¬ 
peak  regimes.  The  choice  of  sign  proves  to  be  very  important  in  calculations,  and  guidelines  for 
sign  change  are  discussed  by  Crisfield  (5).  The  basic  idea  is  that  the  sign  of  the  load  increment 
follows  that  of  the  previous  increment  unless  the  determinant  of  the  tangent  stiffness  matrix 

T 

changes  sign,  in  which  case  a  sign  reversal  is  applied  to  the  load  increment.  Usually,  the  L-D-L 

decomposition  is  used  to  factor  the  tangent  stiffness  matrix.  The  determinant  of  K°  is  equal  to  the 
determinant  of  the  diagonal  matrix,  D,  and  the  appearance  of  negative  eigenvalues  corresponds  to 
the  appearance  of  negative  pivots  in  the  decomposition.  Since  more  than  one  negative  pivot  may 
appear  simultaneously,  the  sign  of  each  pivot,  and  not  the  product  of  the  pivots,  must  be  monitored 
to  detea  the  appearance  of  negative  eigenvalues.  Unfortunately,  numerical  calculations  involve 
finite  digit  arithmetic,  and  if  the  pivots  are  small,  slight  perturbations  may  produce  large  numerical 
errors  (14).  Therefore,  an  incorrect  conclusion  may  be  inferred  with  respect  to  the  sign  of  the  load 
increment,  and  the  iterative  procedure  may  not  converge. 
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To  circumvent  some  of  these  difficulties,  and  to  accommodate  the  physical  feature  that  softening  is 
accompanied  by  localization,  an  alternative  constraint  approach  is  proposed  that  automatically  yields 
the  appropriate  load  increment  With  the  use  of  a  suitable  measure  of  strain,  the  point  in  the  body 
with  the  maximum  value  of  this  measure  is  located  at  each  load  step,  and  an  increment  in  a 
combination  of  the  strain  components  at  this  point  is  prescribed.  The  resulting  equations  represent 
a  slight  extension  to  those  used  previously  by  Crisfield,  de  Borst,  and  their  co-workers  (5, 10). 
However,  as  illustrated  through  the  solution  to  model  problems  in  plane  strain  and  plane  stress, 
such  a  modification  appears  to  result  in  a  solution  scheme  from  which  the  complete  post-peak 
regime  is  easily  obtained. 
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2.0 


THEORETICAL  FORMULATION 


Consider  a  spatially  discretized  system  based  on  the  finite  element  method  in  which  the  vector  of 
nodal  displacements  is  denoted  by  u  and  an  increment  by  Au  .  With  a  superscript  defining  the 
element,  the  incremental  strain-nodal  displacement  relation  over  an  element  is  given  by 


Aee  =  BeAue 


(2-1) 


in  which  e  denotes  the  strain  vector  and  B  relates  the  strain  and  nodal  displacement  components 
through  a  differential  operator  and  interpolation  functions.  If  s  denotes  the  stress  vector  and  C  the 
continuum  tangent  stiffness,  then  the  incremental  relation  between  stress  and  strain  is 


As  =  C:Ae 


(2-2) 


If  R  represents  the  domain  of  an  element,  then  the  element  tangent  stiffness  matrix  is 


Kc=  jBeTCBcdV  (2.3) 

Re 

in  which  dV  denotes  the  volume  element.  The  global  tangent  stiffness  matrix,  K,  is  obtained 
through  a  sum  of  the  element  contributions. 

Only  proportional  loading  will  be  considered;  then  the  load  vector  is  given  by  pq  in  which  q  is 
fixed  and  p  represents  the  magnitude  of  the  load.  For  nonlinear  problems  K  depends  on  u  from 
geometrical  and  material  nonlinearities  of  which  only  the  latter  will  be  considered.  One  method  for 
obtaining  solutions  is  to  solve  the  incremental  equilibrium  relation 

KAu  =  Apq  (2-4) 

in  which  Ap  is  given  and 

p  =  pp  +  Ap  u=up  +  Au  K  =  K(u)  (2-5) 
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The  superscript  p  denotes  the  variable  at  the  end  of  the  previous  incremental  step.  Because  the 
updated  stiffness  matrix  is  not  known,  an  iterative  procedure  is  required  to  obtain  a  solution.  If  I  is 

used  as  the  iteration  count,  and  with  u  0  =  up ,  one  approach  is  to  solve  iteratively 

Ki-i6ul  =  A^  ui=ui-i+5ui  K  i_i  =  K(u  1-1)  (2-6) 

until  a  convergence  criterion  is  met.  However,  for  problems  in  which  a  limit  load  and  softening  are 
expected,  Ap  cannot  be  prescribed  and,  instead,  must  be  determined  indirectly  as  part  of  the 
iterative  solution  procedure.  To  remove  the  resulting  indeterminacy,  a  combination  of  displacement 
components  and  load  increment  is  often  prescribed  in  an  "arc  control"  method,  which  permits  a 
solution  to  be  obtained  even  if  the  tangent  stiffness  matrix  contains  zero  or  negative  eigenvalues.  If 
this  constraint  is  left  open  for  the  moment,  such  procedures  are  obtained  by  recasting  (2-6)  in  the 
form 


K  {  ,  =  5p,q  +  r,_j  r,_,  =  p^q  -  JVs  MdV 

R 

p,  =  p,_i+8p,  p0  =  pp 


(2-7) 


in  which  the  residual  rj_  j  is  included  to  prevent  incremental  drifting  away  from  the  equilibrium 

state.  This  equation  is  identical  to  that  derived  by  de  Borst  (10),  but  with  a  different  point  of  view. 
However,  an  alternative  form  of  (2-7)  is  more  common.  Let 


Pi  =  ^o<1- J®1*  I-ldV 

R 


(2-8) 


Then  (2-7)  becomes 


Ki-i6ui  =  aM  +  Pi-i  P|  =  Po  +  aPi 


(2-9) 


At  this  stage  the  subincrement  in  the  solution  vector  is  considered  to  be  the  sum  of  two  parts: 
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5u  j  =  A^,5u^'  +  8u  [ 


Au0  =  ° 


(2-10) 


Au  j  =  Au  j  j  +  5u , 


K,_,8“l  =  1 


Ki-isur=  P.-1 


With  this  formulation  the  proposed  method  involving  total  strain  control  and  its  relationship  to 
existing  methods  can  be  conveniently  described.  Conventional  arc  control  methods  specify  the  load 
increment  in  terms  of  global  norms  of  displacement  increments  at  the  current  iteration.  Modified 
methods  include  contributions  to  the  arc  control  of  the  increment  in  the  displacement  from  the 
previous  iteration.  The  basic  point,  however,  is  that  global  norms  are  used.  For  problems 
involving  softening  and  localization,  such  norms  may  not  be  sensitive  enough.  Instead  it  is 
proposed  here  that  a  "localized"  norm  be  used  for  this  particular  class  of  problem. 

With  the  separation  of  the  displacement  increments  into  the  parts  given  by  (2-10),  it  follows  from 
(2-1)  that  the  increment  in  the  strain  vector  for  an  element  is 


Aee=BeAup+  ApBcAuq  (2-1  li 

Suppose  a  linear  combination  of  the  increments  in  strain,  Ae,  is  prescribed.  Define  an  element 
constraint  vector  ce  such  that 


Ae  =  cc  Ae' 


(2-12) 


Since  Au  j  =  5u  pf  =  0  for  I  =  1  if  the  convergence  criterion  is  satisfied  at  the  end  of  the  previous 

increment,  and  since  (2-12)  is  to  be  fixed  throughout  the  load  step,  then  the  load  increment.  A  p  { . 
must  satisfy  the  following  equation: 
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Ae 


for  I  =  1 


AjXj  = 


ccV8u« 


AH,  =  - 


ceV8u[ 

c'VSuJ 


for  I  >  1 


(2-13) 


For  a  general  problem,  it  is  logical  to  choose  the  governing  element  to  which  the  constraint  is  to  be 
applied  as  that  element  which  has  experienced  the  most  distortion  as  reflected  through  a  suitable 
norm  of  strain.  For  problems  in  which  softening  is  a  possibility,  the  critical  element  will  be  the  one 
that  first  reaches  the  limit  state,  and  this  element  is  usually  the  one  that  continues  to  be  critical 
because  it  softens  the  most.  The  use  of  (2-13),  with  the  increment,  Ae,  always  assumed  to  be  of 
the  same  sign  for  consecutive  load  steps,  avoids  the  difficulties  associated  with  a  singular  tangent 
stiffness  matrix  at  the  limit  load  and  with  multiple  negative  eigenvalues  in  the  post-peak  regime. 
Rots  (15)  proposed  a  similar  approach  with  a  constraint  condition  involving  the  displacement 

components  on  either  side  of  an  active  crack.  However,  ue  prescribes  Apj  instead  of  the  crack 
opening  displacement  in  the  first  iteration.  As  a  result  there  is  a  potential  problem  involving  the 
monitoring  of  negative  pivots  for  transitioning  critical  points.  With  the  proposed  approach 

involving  total  strain  control,  the  sign  and  magnitude  of  Ap^  are  automatically  adjusted  throughout 
the  iteration  procedure. 

With  regard  to  a  convergence  criterion,  the  results  of  numerical  investigations  indicate  that  the 
energy  criterion,  which  involves  the  product  of  the  out-of-balance  force  with  the  increment  in 
displacement,  is  not  satisfactory,  although  a  component  of  force  may  be  large,  the  corresponding 
component  of  displacement  increment  is  sufficiently  small  to  cause  the  criterion  to  be  satisfied. 

This  situation  arises  frequently  at  limit  points.  Therefore,  in  the  numerical  examples  of  the  next 
section,  convergence  is  based  on  the  following  combined  force  and  displacement  criterion: 
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(2-14) 


where  the  2-norm  is  inferred.  This  criterion  is  also  used  by  Padovan  and  Tovichakchaikul  (9). 
Although  different  values  for  the  tolerances  can  be  specified,  for  this  study  they  were  chosen  to  be 

e  u  =  e  ^  =  0.001 .  For  most  cases  the  part  based  on  displacement  is  satisfied  first. 
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3.0 


NUMERICAL  SOLUTIONS  TO  PLANE  PROBLEMS 


To  illustrate  the  robustness  of  the  proposed  solution  scheme,  a  simple  bilinear  hardening  and 
softening  model  is  used  to  simulate  material  behavior.  The  nonlocal  aspect  of  post-peak  response 
is  not  addressed  and,  therefore,  post-peak  solutions  that  depend  on  mesh  size  are  to  be  expected. 

The  invariants  of  stress,  s  ,  and  plastic  strain,  eP ,  are  defined  in  the  usual  manner  from  the 
deviators  of  stress,  sd  and  plastic  strain,  epd : 

s  =  ytr(sd  •  sd)  ep  =  J|(trdepddepd)  (3-D 

The  yield  function  is  assumed  to  be 


in  which  H  is  a  hardening-softening  function  and  Hl  denotes  the  limit,  or  maximum,  value.  With 
the  definition 


e*= 


(3-3) 


a  bilinear  relation  for  the  hardening-softening  function  is  conveniently  described  as  follows: 


H  =  H0+(HL-H0)fF 

0£ep< 

eL 

H  =  HL(1  -  me*) 

0  £e*< 

x 

ii 

o 

e**_L 

c  *  m 
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_p 

in  which  Hq  denotes  the  elastic  limit.  The  value  of  the  plastic  strain  invariant  at  the  limit  state,  e  L , 

is  assumed  to  be  a  unique  value  no  matter  which  stress  path  is  followed.  The  parameter,  m,  merely 
reflects  the  rate  of  softening.  The  yield  and  consistency  conditions  are  given  by  f  =  0  and  df  =  0, 
respectively.  An  associated  flow  rule  is  used.  Normally,  solutions  to  model  problems  in  plane 
strain  are  given  because  the  displacement  formulation  implies  that  die  plane  strain  condition  is 
satisfied  easily.  However,  with  the  exception  of  some  additional  algebra  involving  the  use  of  the 
tangent  moduli,  solutions  to  problems  in  plane  stress  can  also  be  obtained  for  the  same  geometric 
configurations  with  little  additional  computational  effort.  Both  classes  of  plane  problems  are 
considered  next. 

The  geometry  and  notation  for  the  plane  problem  is  shown  in  Fig.  1.  The  dimensions  used  for  the 
analysis  are  Lx  =  1  m  and  Ly  =  1.5  m.  In  order  to  accommodate  isochoric  deformation,  the 

argument  of  Nagtegaal  et  al.  (16)  for  the  use  of  a  four-node  element  composed  of  four  triangles 
defined  by  the  two  diagonals  in  the  quadrilateral  is  considered  to  be  persuasive  and,  therefore,  the 
same  element  is  used  in  this  study.  The  lateral  surfaces  are  stress-free  (except  for  plane  strain) 
while  the  boundary  condition  at  y  =  0  consists  of  zero  vertical  displacement  and  a  lateral  constraint 
of  zero  displacement  at  one  point  which  is  x  =  0  unless  specified  otherwise.  The  load  is  applied 
through  a  mechanism  that  provides  no  lateral  constraint  and  ensures  that  the  displacement  in  the  y- 
direction  is  the  same  for  all  points  on  the  upper  surface.  Since  the  plasticity  model  is  that  of  von 
Mises,  significant  plastic  deformations  are  exhibited  in  shear  bands  oriented  at  45  deg  to  either 
coordinate.  Another  implication  that  follows  from  the  use  of  the  von  Mises  model  is  that  the  same 
results  are  obtained  for  loading  in  tension  and  compression  so  such  a  differentiation  is  not  made  in 
the  subsequent  description.  Element  mesh  configurations  I,  II  and  in  are  defined  tobe3x3,5x7 
and  9x13  rectangular  grids,  respectively. 

In  the  elastic  regime,  the  material  is  assumed  to  be  isotropic  as  defined  through  Young's  modulus. 
E,  and  Poisson's  ratio,  v .  Material  parameters,  which  are  considered  to  be  representative  of  a  class 
of  geological  materials,  are  chosen  as  follows  (unless  specified  otherwise):  E  =  50  GPa,  v  =  0.2. 

_  p 

e  L  -  0.01 1,  Hq  =  0,  Hl  =  55  MPa,  and  m  =  0. 1 .  The  zero  value  for  the  elastic  limit  is  chosen 

partly  for  convenience  and  partly  so  that  unloading  can  be  easily  detected  because  the  loading  and 
unloading  moduli  are  different  right  from  the  initiation  of  loading.  Initial  imperfections  are 
_P 

introduced  by  letting  e  L  =  0.01  and  H^  =  50  MPa  at  a  point  in  the  body  designated  with  the 
symbol  V. 
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Force 
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The  first  model  problem  is  one  of  uniform  loading  in  the  y -direction  under  the  assumption  of  plane 
strain  (components  with  a  subscript  z  are  zero)  with  all  components  of  displacement  prescribed  to 
be  zero  for  the  point  at  the  origin.  The  constraint  condition  for  the  total  strain  increment  in  the 
critical  element  was  chosen  such  that  the  increment  in  y-y  component  of  normal  strain  is  prescribed. 
The  evolution  of  the  deformation  field  for  element  configuration  HI  is  shown  in  Figs.  2(a)  and  2(b) 
for  plane  strain.  Because  of  the  symmetrical  location  of  the  imperfections  obtained  by  using  an  odd 
number  of  elements  in  each  coordinate  direction,  two  shear  bands  develop  and  intersect  one  end  of 
the  body.  For  the  same  mesh  configuration  and  also  for  plane  strain,  Figs.  3, 4 , 5,  and  6  show 
the  failure  modes  that  develop  for  different  choices  of  the  locations  of  initial  imperfections  and 
slightly  different  boundary  conditions.  The  sensitivity  of  failure  modes  to  imperfections  and 
boundary  conditions  has  been  pointed  out  previously  by  Leroy  et  al.  (17).  For  Figs.  3  and  4,  the 
lateral  constraint  is  at  x  =  0,  while  for  Figs.  5  and  6  the  constraint  is  at  x  =  1.0  and  x  =  0.5, 
respectively.  Figure  3  shows  crossed  shear  bands  while  Fig.  4  shows  a  single  shear  band  that 
does  not  intersect  either  the  top  or  bottom  surface.  Figure  5  is  a  mirror  image  of  Fig.  4  because  of 
the  switch  in  location  of  the  lateral  constraint  on  the  boundary.  The  symmetrical  location  of  the 
lateral  boundary  constraint  in  Fig.  6  is  reflected  in  a  correspondingly  symmetrical  deformation 
field. 

For  plane  stress  and  for  the  point  at  the  origin  fixed.  Fig.  7  illustrates  the  corresponding  final  post¬ 
peak  deformation  field. 

The  effect  of  the  size  of  the  prescribed  total  strain  increment  for  mesh  configuration  I  is  shown  in 
Fig.  8  for  plane  strain.  For  the  largest  increment  of  0.01  the  peak  load  is  achieved  in  only  two 
steps  so  larger  increments  would  not  be  meaningful.  The  fact  that  no  signicant  differences  are 
displayed  can  probably  be  attributed  to  the  use  of  a  bilinear  hardening-softening  function. 

For  plane  strain  and  for  mesh  configuration  I,  load-deflection  curves  are  given  in  Fig.  9  for  various 
values  of  the  softening  parameter.  For  m  =  0.01 , 0.1  and  1.0,  the  peak  load  per  unit  thickness  and 
the  corresponding  deflections  are  (58.6  MN/m,  0.0171  m),  (58.3  MN/m,  0.0160  m)  and 
(52.8  MN/m,  .0140  m),  respectively.  The  loads  and  deflections  shown  in  the  figure  are 
normalized  with  respect  to  these  values. 

For  plane  stress  and  also  for  mesh  configuration  I,  Fig.  10  shows  the  effect  of  the  softening 
parameter  m  on  the  load-deflection  response.  For  m  =  0.01, 0.1,  and  1.0,  the  peak  loads  per  unit 
thickness  and  the  corresponding  deflections  are  (56.1  MN/m,  0.0169  m),  (55.7  MN/m,  0.0168  m) 
and  (53.6  MN/m.  0.0150  m),  respectively.  As  can  be  seen  for  both  plane  strain  (Fig.  9)  and  plane 
stress  (Fig.  10),  numerical  difficulties  are  not  encountered  for  the  full  range  of  post-peak  behavior 
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Deformed  mesh  in  the  post-peak  regime  for  plane  strain  and  constraint  at  x  =  0. 
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Figure  3.  Final  post-peak  deformation  field  with  two  weak  elements  in  the  uppermost  comers. 
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Figure  4.  Final  post-peak  deformation  field  with  the  weak  element  at  the  center. 
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Figure  5.  Final  post-peak  deformation  field  with  constraint  at  x  =  1 . 
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Figure  6.  Final  post-peak  deformation  field  with  constraint  at  x  =  0.5. 


Figure  7.  Final  post-peak  deformation  field  for  plane  stress. 
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Figure  8.  The  effect  of  the  size  of  prescribed  total  strain  increments  on  numerical  solutions 
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Figure  10.  The  effect  of  the  softening  parameter  on  load-displacement  curves  for  plane  stress. 
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ranging  from  almost  perfect  plasticity  to  a  slight  reversal  where  the  post-peak  response  is  taken 
down  to  a  load  that  is  close  to  zero.  However,  it  must  be  emphasized  that  for  this  study  nonlinear 
geometrical  effects  are  not  included,  so  that  the  potential  problems  discussed  by  de  Borst  (18), 
which  are  associated  with  mesh  distortion  and  lockup,  are  precluded. 

T 

In  monitoring  the  negative  pivots  in  the  L-D-L  decomposition  of  the  tangent  stiffness  matrix,  it 
was  noted  that  for  the  case  of  plane  stress,  negative  pivots  did  not  appear  until  those  elements  in  the 
softening  zone  attained  a  high  degree  of  softening.  In  other  words,  it  seemed  that  the  appearance  of 
negative  eigenvalues  (immediate  post-peak  regime)  was  not  accompanied  by  the  appearance  of 
negative  pivots.  Since  this  result  is  contrary  to  basic  theory  in  matrix  analysis,  there  must  be  an 
explanation  involving  round  off  and  numerical  accuracy.  The  results  for  plane  strain  are  consistent 
with  what  might  be  expected  from  theory. 

Because  a  local  constitutive  model  is  used,  a  mesh  dependence  is  expected  for  the  post-peak 
response.  Such  a  dependence  is  shown  in  Figs.  1 1  and  12  for  plane  strain  and  plane  stress  with 
m  =  0. 1  and  a  prescribed  strain  increment  of  of  0.001 .  The  need  for  nonlocal  models  to  prevent 
this  mesh  dependency  has  been  discussed  by  Schreyer  and  Chen  (19),  Chen  and  Schreyer  (20), 
and  Pijaudier-Cabot  and  Bazant  (21),  among  others.  Chen  (22)  has  extended  the  above  analysis 
procedure  to  show  that  the  mesh  dependency  is  significantly  reduced  for  soil-structure  interaction 
problems  when  a  nonlocal  constitutive  model  is  used. 


22 


Figure  12.  Illustration  of  mesh  dependence  for  plane  stress. 
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4.0  CONCLUSIONS 


A  slight  modification  to  existing  arc-control  methods  for  obtaining  numerical  solutions  to  problems 
exhibiting  post-peak  softening  is  proposed  in  which  a  particular  linear  combination  of  the 
components  of  the  total  strain  increment  is  prescribed  at  a  critical  point.  There  is  no  reason  why  the 
method  could  not  be  easily  extended  to  include  any  norm  of  the  total  strain  increment  tensor.  Then 
the  load  increment  is  automatically  determined  with  the  result  that  the  procedure  is  particularly 
robust  for  transitioning  critical  points  at  which  an  eigenvalue  of  the  tangent  stiffness  matrix  is  zero. 
The  examples  involve  materially  nonlinear  problems,  but  the  procedure  should  also  be  applicable  to 
geometrically  nonlinear  cases. 

The  examples  show  that  the  evolution  of  failure  modes  is  quite  sensitive  to  the  location  of  initial 
imperfections  and  to  boundary  conditions.  Such  results  are  in  agreement  with  many  experiments 
involving  materials  such  as  concrete  where  laboratories  have  reported  significantly  different  failure 
loads  and  failure  mechanisms  for  specimens  composed  of  one  batch  of  concrete.  The  current 
analysis  indicates  that  such  diversity  may  be  entirely  consistent  with  the  slightly  different  boundary 
conditions  associated  with  different  experimental  apparatuses. 
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DISTRIBUTION  LIST 


Dept  of  the  Air  Force 


AUL/LSE  1  copy 

Maxwell  AFB ,  AL  36112 


KQ  AU 

Attn:  Tech  Library 
Maxwell  AFB,  AL  36112-5001 


AF  Institute  of  Technology/DE 
Attn:  Tech  Library 
Wright; Patterson  AFB,  OH  45433-5001 


FCDNA/FC 

Kirtland  AFB,  NM  87115-5000 


1  copy 


1  copy 


1  copy 


USAFA  1  copy 

Attn:  Tech  Library 

Dept  of  Civil  Engineering 

Colorado  Springs,  CO  80840-5001 


BMO/MYET 
Norton  AFB, 

CA 

92409-6468 

1 

copy 

HQ  AFESC 

Attn:  Tech  Library 

Tyndall  AFB,  FL  32402-6001 

1 

copy 

HQ  AFESC/RD 
Tyndall  AFB, 

FL 

32542-5000 

1 

copy 

HQ  AFSC/XTH 
Andrews  AFB, 

MD 

20334-5000 

1 

copy 

HQ  AFSC/DE 
Andrews  AFB, 

MD 

20334-5000 

1 

copy 

HQ  PACAF/DE 
Hickara  AFB, 

HI 

90853-5001 

1 

copy 

HQ  USAF/LEED 

1 

/ 

copy 

foiling  AFB 

Washington  DC  20330-5120 

1  copy 


Weapons  Laboratory 
Attn:  NT 

Kirtland  AFB,  NM  87117-6008 
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DISTRIBUTION  LIST  (CONTINUED) 


Weapons  Laboratory 
Attn:  NTE 

Kirtland  AFB ,  NM  87117-6008 

Weapons  Laboratory 
Attn:  NTED 

Kirtland  AFB,  NM  87117-6008 

4  Weapons  Laboratory 

Attn:  NTEO 

Kirtland  AFB,  NM  87117-6008 

Weapons  Laboratory 
Attn:  NTES 

Kirtland  AFB,  NM  87117-6008 

Weapons  Laboratory 
Attn:  NTESG 

Kirtland  AFB,  NM  87117-6008 

Official  Record  Copy 
WL/NTESR  (Maj  John  Gill) 
Kirtland  AFB,  NM  87117-6008 

Weapons  Laboratory 
Attn:  SUL 

Kirtland  AFB,  NM  87117-6008 

Weapons  Laboratory 
Attn:  HO 

Kirtland  AFB,  NM  87117-6008 


Dept  of  the 


Harry  Diamond  Lab 
Attn:  Library 
2800  Powder  Mill  Road 
Adelphi ,  MD  20783 

US  Army  Engr  Waterways  Exper  Sta 
Attn:  Tech  Library 
P.0.  Box  631 
Vicksburg,  MS  39181 


1  copy 

1  copy 

1  copy 

2  copies 

2  copies 

2  copies 

2  copies 

2  copies 

Army 
1  copy 

1  copy 
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DISTRIBUTION  LIST  (CONTINUED) 


US  Army  Engr  Waterways  Exper  Sta  1  copy 

Attn:  Dr  Rohani 
P.0.  Box  631 
Vicksburg,  MS  39181 

US  Army  Construction  Engineering  1  copy 

Research  Laboratory 
Attn:  Technical  Library 

P.O.  Box  4005  %■ 

Champaign,  IL  61820-1305 

US  Army  Corps  of  Engineers  1  copy 

Omaha  District 

MORED-ST 

6014  USPO  &  Courthouse 
Omaha,  NE  68102-4978 


Dept  of  the  Navy 


Naval  Research  Lab  (NRL)  1  copy 

Attn:  Code  2627,  Tech  Lib 

Washington  DC  20375 

Naval  Surface  Weapons  Center  1  copy 

Attn:  Tech  Library 
Dahlgren,  VA  22448 

Naval  Civil  Engrg  Lab  1  copy 

Attn:  Library  (Code  L08A) 

Port  Hueneme ,  CA  92043 

Naval  Surface  Weapons  Center  1  copy 

White  Oak  Laboratory 

Attn:  Code  X221,  Tech  Library 

Silver  Spring,  MD  20910 


Government  Agencies 


Defense  Nuclear  Agency  1  copy 

Attn:  SPSD  (Kent  Goaring) 

Washington  DC  20301 

DTIC/FDAC  2  copies 

Cameron  Station 
Alexandria,  VA  22304-6145 
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DISTRIBUTION  LIST  (CONCLUDED 


Director 

Attn:  SPSD  (LTC  Rubio) 
Defense  Nuclear  Agency 
Washington  DC  20305-1000 


Government 


AFELM 

Rand  Corporation,  Library  D 

P.O.  Box  2138 

Santa  Monica,  CA  90406 

Sandia  National  Laboratories 
Attn:  Tech  Lib,  Div  3141 
Kirtland  AFB ,  NM  87185 

University  of  New  Mexico 
Attn:  D.  Calhoun 

NMERI 

P.O.  Box  25,  University  Station 
Albuquerque,  NM  87131 

Karagozian  &  Case 
Attn:  John  Karagozian 

620  N.  Brand  Blvd  #300 
Glendale,  CA  91203 


1  copy 


Contractors 


1  copy 


1  copy 


1  copy 


1  copy 
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